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ABSTRACT 

We introduce and implement two novel ideas for modeling lensed quasars. The 
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I first idea is to require different lenses to agree about Hq. This means that some 

models for one lens can be ruled out by data on a different lens. We explain using 
^ I two worked examples. One example models 1115+080, 1608+656 (time-delay quads) 

' and 1933+503 (a prospective time-delay system) all together, yielding time-delay pre- 

^ ■ dictions for the third lens and a 90%-confidence estimate of -f^o"^ = 14.6l^'y Gyr 

{Ho = 67 +26 km s-^Mpc-i) assuming {Qm = 0.3, = 0.7). The other example 
models the time-delay doubles 1520+530, 1600+434, 1830-211, and 2149-275, which 
gives -f^o"^ = 14.5+^'5 Gyr {Hq = 67+^^3 km s^^Mpc^^). Our second idea is to write 
the whole modeling software as a highly interactive Java applet, which can be used both 
for coarse-grained results inside a browser and for fine-grained results on a workstation. 
Several obstacles come up in trying to implement a numerically- intensive method thus, 
but we overcome them. 



Subject headings: gravitational lensing 
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1. Introduction 

Some aspects of modeling lensed quasars are much as they were just after the discovery of 
the first double quasar 0957+561. In one of the earliest lens-modeling papers, Young et al. (1981) 
are concerned with some now- very-familiar issues: the effect on image positions of both the main 
lensing galaxy and other galaxies, the time delays predicted by the models, the non-uniqueness of 
the models despite the adequacy of the data, and the desirability of supplementary data about the 
lens, such as velocity dispersions or X-rays. 

But other aspects of the subject these days would have been unimaginable in 1981. The first 
double quasar has been joined by dozens of others: the CASTLES survey compilation (Kochanek 
et al. 1998) currently lists 76 secure multiple-image galaxy-lens systems, with image positions 
measured to the mas level, even more precisely if there are compact radio sources. And the time 
delay for 0957-1-561, for which Young et al.'s preliminary estimate was about 5 years, is now 
measured as 423 it Idays (Oscoz et al. 2001), along with time-delay measurements for eight other 
systems (Schcchter et al. 1997; Lovell et al. 1998; Biggs et al. 1999; Burud et al. 2000; Cohen et al. 
2000; Burud et al. 2002a,b; Fassnacht et al. 2002; Hjorth et al. 2002). 

These excellent data demand new, more automatic, and more portable software tools for 
modeling the lenses. Thus motivated, we have developed a new code, PixeLens, which we present 
in this paper. 

PixeLens works by reconstructing a pixclated mass map for the lens, an idea we first im- 
plemented in Saha &: Williams (1997). Most lens modeling codes work by fitting a parametric 
functional form — for example. Young et al. (1981) fitted King models; gravlens by Keeton (2001) 
is a modern example, offering the user a large choice of parametric models. There are other possi- 
bilities too: Trotter et al. (2000) reconstruct lenses non-parametrically as we do, but use multipole 
expansions rather than pixels. 

PixeLens generates large ensembles of models rather than one or a few mass maps, as a way of 
addressing the non-uniqueness problem. We introduced this strategy in Williams h Saha (2000) for 
pixelated models. Keeton &; Winn (2003) use a somewhat similar strategy, but with parametrized 
models. 

But PixeLens also brings two completely new features, one astrophysical, one computational: 

• It can model several lenses simultaneously, enforcing consistency of Hq across different time- 
delay lenses. As a result, lenses can be used to constrain other lenses in an interesting and 
surprising way; 

• The code is highly portable and can run without change as a standalone program, or as a 
Java applet inside a web browser. In the online version of this paper, it is available as an 
alternative version of Figure 1. 
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The best way to explain what we have implemented and what problems remain is through an 

example. So in the following section we will work through a simultaneous reconstruction of three 
lenses: the time delay quads 1115+080 and 1608+656, and the ten-image system 1933+503. This 
will be the main part of the paper. We turn to some computational issues after that. 

2. A worked example: reconstructing three lenses 

In our first worked example, we consider the three lenses currently best-constrained by obser- 
vations. 

1. 1115+080 is a quad where the lens is an elliptical galaxy supplemented by external shear 
from other group galaxies. It was the second lens to be discovered (Weymann et al. 1980). 
Schechter et al. (1997), together with an improved time-series analysis by Barkana (1997) 
provide time-delays between two pairs of images. 

2. 1608+656 is a quad with the lens being apparently two interacting galaxies. It was discovered 
(Myers et al. 1995) in the CLASS survey. There are time delays between all pairs of images 
(Fassnacht et al. 2002). 

3. 1933+503 is really two quads and a double. The lens appears to be an elliptical, but the source 
is a radio source with a core and opposing lobes. The core and one of the lobes are imaged 
as quads while the other lobe is imaged as a double. It was also discovered in the CLASS 
survey (Sykes et al. 1998). There are no measured time-delays, but there are prospects for 
time delays from the core quad. We will predict time delays for the core quad. 

The CASTLES web page provides further information. We can also make some preliminary infer- 
ences about the lenses by examining the morphology, as explained in Saha &: Williams (2003): we 
can work out the time-ordering of each of the image systems, and for 1115+080 we can recognize 
an external shear and its approximate axis. 

Let us now proceed to apply PixeLens. Figure 1 shows the user interface, which we will explain 
in stages below. 

2.1. The data input 

Data is entered by typing into the 'data area'. (There is also an 'input paste' feature with a 
number of example inputs, which can be pasted and then edited by selecting the 'edit' option of 
the menu.) 

For 1115+080, we enter the following. 
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object 1115 

symm pixrad 12 maprad 2 
redshifts 0.311 1.722 
shear -45 
quad 

0.3550 1.3220 
-0.9090 -0.7140 13.3 
-1.0930 -0.2600 

0.7170 -0.6270 11.7 

The first line just labels the object. The second line specifies that the mass maps be inversion- 
symmetric with a radius of 2" and 12 pixels (thus setting the pixel size). The fourth line gives 
the approximate shear axis ('shear 135' would be equivalent), in degrees relative to East. The 
program will consider positions angles within 45° of what is given. The last four lines give the image 
positions (in arcseconds and relative to the lens center) and the time delays (in days and relative 
to the previous image). The images must be ordered by arrival-time. The '0' in the second-last 
line means that the time-delay is > but unknown. Otherwise, uncertainties in image positions 
and time delays are assumed negligible. 

Note how the input format assumes that some preliminary analysis, of the qualitative variety 
explained in Saha k, Williams (2003) has already been done. 

For 1608-1-656 we enter the following. 

object 1608 
pixrad 9 maprad 2 
redshifts 0.630 1.394 

quad 

-1.300 -0.800 
-0.560 1.160 31 
-1.310 0.700 5 
0.570 -0.080 40 

The format is just as before, except that we have left out symm, because the lens here appears to 
be asymmetric. 
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For 1933+503 we enter the following. 

object 1933 
symm pixrad 10 
redshifts 0.76 2.63 
shear 45 

quad -0.40 0.53 0.43 -0.26 

0.40 0.19 -0.19 -0.33 
quad 0.54 -0.44 -0.15 0.44 
-0.03 0.45 -0.38 -0.12 
double -0.47 0.60 0.13 -0.30 

Here we have left out the mass-map radius, leaving the program to choose it. We have also given 
image data on two quads and a double rather than just one quad, and indicated unknown time 
delays with many Os. Line breaks in the input have no significance. 

Figure 2 summarizes the image data and the pixelation. (For symmetric mass maps, only half 
the pixels are shown here, the other half are specified by inversion symmetry.) PixeLens produces 
such plots before starting the real computations, and they provide a convenient goof test on the 
input. 

Now we type 
models 100 

in the 'data area' and click on the 'run' button, and away we go. 



2.2. Generating an ensemble of models 

At this point, we recall some lensing theory. For a source at (3 , the arrival-time at 9 is 




where k is the dimcnsionless density (also called convergence). If we measure angles in arcsec, r 
will have units of arcsec^. To turn it into days we have to multiply by gT{zi^, za,), where T is a 
time scale and g is Hq^ in Gyr (see Appendix A). In addition to the indicated redshift dependence, 
T{zi,,zs) has a weak cosmology dependence. By default PixeLens assumes JIm = 0.3, J^a = 0.7, 
which gives T^ii^ = 3.34 days, T^^os = 10.28 days, and T^^^s = 9.72 days. 

In our models, the mass is pixelated. Therefore we can simplify Equation (1) to 

T{e) = 0\''-e-p-Y,'^nQnio) (2) 

n 
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where is the density on the n-th pixel and Qn^S) is the integral from (1) evaluated over that 
pixel. The explicit form of QniO ) is a messy but straightforward mixture of logs and arctans and 
is given in Saha &: Williams (1997). The important thing about Equation (2) is that it is linear in 
the unknowns, /3 and Kn- We also allow a constant external shear 

Mol-ol) + '^i20.ey (3) 

to be added to Equation (2); again, the shear is linear in the unknowns 71,72- 
We now implement the following observational constraints. 

1. An image observed at 9i implies 

Vr(^,) = 0. (4) 

In a multiple-image system all the images have the same unknown (3 , so each multiple- 
image system introduces 2 ([images] — [sources]) constraints. Thus a double provides 2 image 
constraints, a quad 6, while 1933-1-503 provides 14. 



2. A time-delay measurement between images at 9i and 9j implies 

[obs delay] 

Here g is another unknown and, fortunately for us, g^^ appears linearly. If a time-delay is 
unmeasured, but we are sure about the ordering, we replace (5) by an inequality constaint 

TiOi) - T(9j) > 0. (6) 

We have 3 equality constraints from time delays in 1608-1-656, 2 in 1115-1-080, and none in 
1933-F503. 

3. We require the images to have the expected parities, and for the image elongation when 
projected along the radial direction to be between and 10. Although very modest, these 
requirements do suppress a tendency for spurious images to be created in the vicinity of a 
nearly-merging pairs of genuine images. Constraints of this type are also linear (AbdelSalam 
et al. 1998). 

4. We remarked after Equation (5) that g is an unknown variable in the constraint equations. 
Now, if we write the constraint equations for several lenses, then each lens will introduce a 
separate unknown variable for g, say g^^^^, g^^'^^ and so on. But g is universal, so we have 
additional constraint equations of the type 

^1115^^1608^ ^1608 ^^1933 (7) 

Though g is still unknown, it now couples the constraint equations for different lenses. This 
coupling is the central astrophysical contribution of this paper. In particular, it means that 
if one lens rules out a certain range of g, PixeLens will not consider models of any lens that 
require g in that range. 
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We also impose secondary constraints based on our general knowledge of galaxy mass profiles. 

1. All the tin > 0, 

2. If the galaxy does not appear very asymmetric, the mass profile is required to have inversion 
symmetry about the lens center. 

3. The density gradient anywhere must point within < 45° of the lens center. The precise 
formulation of this constraint is given in Saha Sz Williams (1997). 

4. The Kn of any pixel must be < 2[average of neighbors], except for the central pixel, which is 
allowed to be arbitrarily dense. 

5. The radial mass profile must be steeper than \6 \~^'^. The stellar dynamical evidence indicates 
that the total density distribution in the central regions of ellipticals is well approximated 
by isothermal, (Rix et al. 1997; Gerhard et al. 2001), while the study of the dynamics 
of the gas in the center of our Galaxy show that total density scales as r~^-'^^ (Binney et al. 
1991). Guided by these studies we choose \9 |~°'^ as a conservative minimal steepness of the 
2D density profiles in galaxies. 

In Bayesian terminology, these secondary constraints constitute a prior. 

The observational and prior constraints confine allowed lenses to a polyhedron in the space 
of K2, . . . , kn, P )• PixeLens searches this polyhedron by a Monte-Carlo technique, briefiy ex- 
plained in Figure 3. The key idea is to random- walk between models such that probability of moving 
from model A to model B equals the probability of moving from model B to model A. This is a 
simple example of a Metropolis algorithm, or a Markov-chain Montc-Carlo method. Generalizations 
to non-uniform priors are also possible — see, for example, Saha (2003). 

In our example, each model that PixeLens generates is a three-lens model: there arc mass 
maps, external-shear models where required, and inferred source positions for each of 1115-1-080, 
1608+656, and 1933+503. In any one model the lenses are coupled by a common value of g, but 
different models may have different values of g. 

As soon as the first model is ready, the user can examine it. But the users will not want to 
examine every model in an ensemble of 100. So PixeLens presents two kinds of results from the 
ensemble: details of the average model of the ensemble-so-far and various statistics derived from 
the ensemble. (By ensemble-average model we mean a model with the densities in each pixel, the 
source positions, and any external shear averaged over the ensemble. The constraint equations 
being linear in the quantities being averaged, they are still satisfied by ensemble averages.) Let us 
consider the two kinds of results in turn. 
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2.3. Results: the ensemble-average model 

Figure 4 shows mass maps for the three lenses from the ensemble-average model. For 1115+080 
and 1608+656, the mass maps are very similar to those in Williams & Saha (2000). The modeling 
method and input data in the earlier paper were very similar — except for the coupling (7) — but the 
code was completely different. We remark that the elongation of the mass profiles in 1115+080 is 
towards the other group galaxies, the elongations in 1608+656 and 1933+503 are what we might 
expect from the light in those lenses. 

Figure 5 shows the lens potential in 1115+080. (The default is to pick the contours that pass 
through the images, but the user can try different contour spacings interactively, by typing in the 
appropriate values in the 'plot parameters' boxes.) Notice how smooth the potential is: nothing of 
the pixelation so evident in the mass maps is visible here; it has been washed out by the double 
integral in Equation (1). 

Figure 5 also illustrates a convenient way of gauging the importance of external shear in a 
model. PixeLens allows the user to plot the lens potential with the external shear exaggerated. 
(Such exaggeration is purely for visualization, it has no effect on the models.) The user can adjust 
an exaggeration parameter: the default is 1, while means omit the external shear. Figure 5 shows 
that a 3-fold exaggeration makes a drastic difference to the potential, indicating that the external 
shear in this lens is large but not dominant. 

Figure 6 show saddle-point contours in the arrival-time surfaces for the two quads and a double 
of 1933+503. Saddle-point contours are the default, but as with the potential, the user can try 
different contour spacings interactively. Examining the arrival-time contours is a very good way of 
checking that the model is sensible: small errors in the input or the modeling code nearly always 
lead to noticeable spurious features in the arrival-time surface, such as contorted contours or extra 
images. 

Arrival-time contours also have another, somewhat surprising, use. Figure 7 shows the result 
if the user specifies a contour-spacing of 0.004. (The units here are arcsec^. We can convert to 
physical time units, as explained after Equation 1, and the result is about 5hr, but let us continue 
with arcsec^.) A plot of the arrival-time surface with closely-spaced contours tends to resemble an 
Einstein ring. In fact, it approximates — see Saha & Williams (2001) for details — the image of a 
source with a conical light profile of radius 

[contour spacing] 
[thickness of contour lines] ' 

Now, the line thickness in all the plots in this paper is ^ of the plot width. For Figure 7, that 
amounts to ~ 0.015". Hence Figure 7 predicts the image due to a conical profile of radius ~ 0.33" 
centered on the QSO. In fact it resembles the ring imaged by Impey et al. (1998), though there are 
differences of detail. Thus, we can make quantitative inferences about the host galaxy simply by 
inspecting images and doing some mental arithmetic. 
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2.4. The model ensemble: stats 

Out of the full ensemble of models, it is possible to extract many kinds of statistics. PixeLens 
implements four kinds, all involving Hq^ . 

Figure 8 shows a histogram of g values from the ensemble, and has a straightforward interpre- 
tation as the posterior distribution of H^^. This figure is analogous to Figure 17 in Williams & 
Saha (2000), but better. In the earlier work we had separate model ensembles for 1115+080 and 
1608+656 and multiplied their /t-histograms bin-by-bin. That procedure had the problems of (i) a 

reduced number of models in the common region, and (ii) increased shot noise from bin-by-bin 
multiplication. In PixeLens, because of the constraint (7), there is only one histogram of g, which 
eliminates both problems. This histogram (or rather, the unbinncd values it represents) gives 

Hq-^ = 14.6+i;2 Gyr {Hq = 67+^4 local units) at 68% confidence 
i?o"^ = 14.6+?:^ Gyr {Hq = 67+26 local units) at 90% confidence. 

The median Hq is somewhat higher than in Williams k. Saha (2000), but this is expected 

because that paper used an Einstcin-de-Sitter cosmology, which shortens Ti^zi^^z^) by ~ 5% for 
1115+080 and by ^ 10% for 1608+656. The uncertainties are similar. 

Figure 9 shows the predicted time-delays for the core-quad in 1933+503. 

Figure 10 shows the radial index a from a radial fit of \0\~°^ against g"^ . Most lens systems 
when modeled independently show a correlation between a and g~^ (or /i), in the sense that steeper 
density profiles result in higher values of h (e.g.. Figure 13a and 16a of Williams & Saha 2000). For 
the three coupled lenses shown here, 1115+080 shows the correlation weakly and the others hardly 
at all. 

Figure 11 shows the mean k in the annulus between the innermost and outermost images 
(say («:)ann) against g^^ . Again, when modeled independently, most lens systems show an anti- 
correlation between (K)ann and g^^ , whereas two of the three panels of Figure 11 show no trends 
at all. This can be expected, because coupling different time-delay systems weakens correlations: 
coupling reduces the range of g explored, and over a reduced range of g any correlation involving 
g will appear weaker. For an extreme case, suppose we artificially fix the value of 5, which we can 
do by adding (say) 

g 14 

at the end of the input. Then Figures 10 and 11 will consist of vertical clusters of point at 5 = 14. 

However, the above argument does not fully account for the lack of trends in the middle and 
lower panels of Figure 11; other causes are partially responsible. Kochanek (2002) derives a linear 
relation between time delay, ( 

/i)ann) and /i, as a leading order approximation from lensing theory. If 
one of these three (sets of) variables is fixed, say if the time delays are fixed by observations, then a 
well defined trend would be expected to relate the other two ((K)ann and h). This is demonstrated 
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by 1115+080. 1933+503, on the other hand, is not expected to exhibit such a trend because no time 
delays exist for this system. Finally, 1608+656 shows no trend probably because it is a distinctly 
asymmetric lens. 

2.5. Artifacts 

To finish this worked example, we remark on some artifacts evident in pixelated models. 

The most obvious spurious feature is the pixel-scale substructure. In principle, one could 
replace mass pixels with basis functions, thereby preserving the linearity of the constraint equations 
and leaving the main algorithms unchanged, but making the mass maps smooth. That would be 
a useful development for microlensing applications, where local density and shear are crucial. But 
for studying image positions and time delays, basis functions are unlikely to make any difference, 
because the pixel scale structure is completely suppressed in the lens potential and arrival time — see 
Figure 5. 

Pixelation has a serious side-effect, however, if the pixels are made too large. We have found 
from experimenting with different pixel sizes that if the pixels are too large and the mass-distribution 
is poorly resolved, fewer steep mass-models arc allowed. As a result, the prior is shifted to favor 
lower values of g~^, and hence the Hubble-time histograms also shift towards larger Hq^. The 
effect is most noticeable in the a-g~^ correlation. The pixel-sizes used for this paper are chosen 
small enough not to cause this particular problem. 

A second artefact concerns the central image. Now the central maximum itself is obvious from 
arrival time contour maps. What is further often evident is that the central maximum is slightly 
offset from the lens center and not as steep as expected from the practical unobservability. This 
happens because our mass maps lack a central density cusp, instead distributing the mass over the 
central pixel. In principle, we could incorporate central density cusps into PixeLens. But it would 
make no difference to the observed images because, as is well known, a circular mass disk has the 
same lensing effect outside itself as a point mass. 

A third (and rare) artefact is spurious extra images. The density- gradient constraint, we have 
found, is generally enough to suppress extra images, but it is possible that a few rogue models with 
spurious extra images are present in the ensembles. In principle we could reject those models by 
examining every single arrival-time surface by eye, but we have not implemented that. 

3. A complementary worked example: four lenses 

In our second example we will be more brief. 
We reconstruct four time-delay doubles. 
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• 1520+530 was discovered by Chavushyan et al. (1997) and Burud et al. (2002b) measured its 
time-delay. 

• 1600+434 was discovered by Jackson et al (1995) and Burud et al. (2000) measured its time 
delay. 

• 1830-211 was discovered by Rao & Subrahmanyan (1988) and Lovell et al. (1998) measured 
its time delay. 

• 2149-275 was discovered by Wisotzki et at (1996) and Burud et al. (2002a) measured its time 
delay. 

3.1. The data input 

We have already explained the input format, and without further ado here are the data we 
input. 

object 1520 

symm pixrad 8 redshifts 0.71 1.855 shear 90 
double 1.141 0.395 -0.288 -0.257 130 

object 1600 

symm pixrad 8 redshifts 0.42 1.59 
double 0.610 0.814 -0.110 -0.369 51 
object 1830 

symm pixrad 8 redshifts 0.89 2.51 
double -0.50 0.46 0.15 -0.26 26 

object 2149 

symm pixrad 8 redshifts 0.490 2.03 
double 0.736 -1.161 -0.173 0.284 103 

We can omit the number of models since 100 is the default. 

In the case of PKS 1830-211, the position of the lensing galaxy at z = 0.89 remains uncertain. 
Various scenarios have been proposed in the literature (Lehar et al. 2000; Winn et al. 2002; Courbin 
et a. 2002; Koopmans &: de Bruyn 2003); we use the astrometry from Lehar et al. (2000). 

3.2. Results 

Since readers will have already seen the various kinds of results possible, we now give only a 
small selection. 

Figure 12 shows the ensemble-average mass maps. 
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Figure 12 shows the common histogram of Hubble times, from which we obtain 

Ho~^ = Gyr {Ho = Gltj local units) at 68% confidence 

Ho~^ = U.btli Gyr {Hq = 67+^1 local units) at 90% confidence 

Figure 14 shows ( 

^)ann against g ^ for the four doubles. All four show a correlation, 1600-1-434 
being the best. These correlations are somewhat worse compared to what they would have looked 
like if the lens systems were modeled independently (as the reader can readily verify using PixeLens): 
as explained in Section 2.4, coupling lens systems weakens the correlations. 

4. Computational issues: reconciling portability and efficiency 

The first decision we made about PixeLens was to write it in Java. The overriding consideration 
was portability: we wanted to code to run without modification on a high-end workstation, a laptop 
computer used during a talk, or inside a web browser. While developing the software we also found 
other advantages of Java. For example, the graphics libraries are part of the standard, making 
GUIs very easy to write. Also, the language is remarkably clean^ and easy to debug. 

But these advantages come at a price, and we found there were four issues we had to confront. 

First, Java still lacks compilers optimized for numerical work. Java numerical code runs slower 
than C-|--|- by a factor of 3 or more, and also takes up more memory. This is an annoying in- 
convenience, but it is not fatal. (Given Moore's law, it is like running C-|— |- on a two-year-old 
computer.) 

Second, there is very little support for scientific programming. In particular, in all the elaborate 
Java graphics libraries, there is no simple package to do a plot of x against y\ In response, we just 
wrote whatever utility software we needed. For example, we wrote Java code to generate PostScript 
code to draw the big and small ticks in an x, y plot. But this utility code is now available for other 
projects and to other researchers. 

Third, disk input/output is not allowed from inside a web browser. This is an essential security 
precaution — clearly, untrusted code running over the internet must not be allowed disk access. 

Fortunately, the language allows an elegant workaround. PixeLens checks whether it is running 
inside a web browser or as a standalone program (as an applet or an application in the jargon) and 
activates disk input/output only in the latter case. So users can try the code out over the internet 
inside a web browser; if they are interested enough to want disk output, they can download and 
run like a normal program. 

Fourth, and most importantly, there is no batch mode inside a web browser. If the program 
takes more than a few seconds to produce results, which user will want to wait in front of their web 



^Apart from a culture of capital letters inside words, which made the name PixeLens inevitable. 
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browser? This seems to us the fundamental reason why numericahy-intensive Java programs were 
previously unknown. One can write such programs and run them in batch mode in the normal 
way, but if they arc useless inside a web browser the single largest advantage of Java is lost. But 
there is a way out. That way is to provide the user with intermediate results, and the means to 
post-process them in non-trivial ways, while the main computation is still running. With PixeLens, 
the user can start examining the results as soon as the first model of the ensemble is complete, 
even though only 1% of the total computation may have run. We found that providing the user 
with intermediate results and sufficient flexibility to make them interesting influenced almost every 
aspect of our program design. 

PixeLens is internally documented in a literate programming idiom (Knuth 1992) and we will 
not enter into minutae here, but Figure 15 gives a general idea of how the code is organized. Imagine 
each of the labelled ovals as consisting of many variables and several functions that can operate on 
those variables. Functions from one oval may know something about the variables and functions of 
another oval, but only the bare minimum needed to work together. 

(Text I/O) is the simplest of the ovals. It accepts the data and other input like 'run' and 
'pause', and outputs diagnostics and error messages, but does not interpret the input itself. (Router) 
interprets enough of the input to work out how many lenses there are, allocates the required number 
of (Lens) ovals — two in Figure 15 — and passes the data to them. The (Lens) ovals are the only 
parts that know about lensing theory and pixelation. Each of them turns the data into constraint 
equations. (Router) then takes the constraint equations from each (Lens), packs them into one 
large set of equations, and gives them to (Simplex). (Simplex) generates solutions using linear 
programming and Monte-Carlo, and it does nearly all the numerical work involved. Each time a 
new solution is ready, (Simplex) passes it to (Router), which unpacks it and distributes among the 
(Lens) ovals. Each (Lens) oval then post-processes the solution, preparing grids of k and t{9 ) and 
so on. The post-processed material is passed via (Router) to (Graphics I/O). (Graphics I/O) can 
interact with the user and do some non-trivial computation on its own: it can make contour plots, 
scatter plots, and histograms to the user's specification, all without needing to know any lensing 
theory. Moreover (Graphics I/O) can work without interrupting (Simplex) or the other ovals, even 
for error-recovery. In effect, model generation happens in the background while the user examines 
intermediate results in the foreground. 



5. Conclusions 

This paper continues our work on modeling lens quasars using pixelated mass maps. We elab- 
orate on ideas introduced in earlier papers (Saha &: Williams 1997; Williams & Saha 2000; Saha &: 
Williams 2001; Raychaudhury et al. 2003): (a) formulating lens reconstruction as an undetermined 
linear inverse problem, (b) searching through model space to infer Hq (or alternatively, predict time 
delays) together with systematic uncertainties, (c) predicting the main features of any Einstein ring, 
and if observed, estimating the size of the host galaxy. But more importantly, we introduce and 
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implement two new ideas: (i) reconstructing several time-delay lenses simultaneously while requir- 
ing them to agree about Hq, (ii) making the software so portable that it can be run inside a web 
browser while reading this paper. 

The results in this paper are from two worked examples of seven lenses in all. 

1. First we model the time-delay quads 1115+080 and 1608+656, and the 'dec' 1933+503. The 
two time-delay lenses help check against our earlier work using a different code. By coupling 
them with 1933+503 we can predict time-delays for the latter that are not conditional on Hq. 
We also obtain 

Hq-^ = U.etl:^ Gyr {Ho = 67+26 local units) at 90% confidence. 

2. Then we model the time-delay doubles 1520+530, 1600+434, 1830-211, and 2149-275, and 
here the main result 

Ho~^ = 14.5l^J Gyr {Hq = 67+i| local units) at 90% confidence 

The reference cosmology has {Qm = 0.3, Oa = 0.7). Note that the two worked examples contain 
no lenses in common, so the close similarity of the two derived values of i^o (1 and 2 above) is very 
encouraging. 

The well-known correlations — (i) steeper lenses give higher Hq, and (ii) more mass in the 
annulus covering the images gives lower Hq — are present for most of the lenses, and of these (ii) 
is better. However, such simple correlations appear to get weaker as Hq is better constrained, 
which happens as a consequence of coupling lenses.^ Furthermore, lenses with no time delays, 
like 1933+503, as well as lenses with very asymmetric mass distributions, like 1608+656, are not 
expected to show the above correlations. 

At this stage wc remain somewhat cautious about the Hubblc-timc estimates. As wc explained 
in subsection 2.2, the distribution of models we obtain is conditional on the prior we use. In the 
past, when the uncertainties on the Hubble time were very large, fine-tuning the prior was not 
so important. But now, with the uncertainties shrinking down from better data and improved 
modeling, improving the prior is probably the next priority. 

A. Units 

In lensing the Hubble time Hq^ plays a more fundamental role than the Hubble constant Hq. 
Note that the crucial observable has the same dimensions as H^^ even though its value is many 

orders of magnitude different. 



famous comment by R.O. Redman [quoted by Hcggie & Hut (2003)] comes to mind: hearing "after all, a star 
is a pretty simple thing," Redman retorted "at a distance of 10 parsecs you'd look pretty simple". 
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Accordingly, rather than using the traditional dimensionless number h we adopt a dimension- 
less number g, which is defined by 

H^^=gGyv. (Al) 

Numerically, g = 9.78/h. 

To obtain convenient units for the computations, we define 

T(zl,zs)^(1 + ^l)^^ (A2) 

where the distances are computed using a fiducial value Hq^ = 1 Gyr (and in some chosen cos- 
mological model). Multiplying further by 3.6525 x 10^^/206265^ converts T{zi„zs) to units of 
days arcsec"^. Then, by using T{zi,,zs) in appropriate places we can leave all positions in arcsec 
and all time delays in days. 

We remark that in currently-favored cosmological models, Hq^ is nearly equal to age of the 
universe. In a flat cosmology containing only matter and dark energy, the Hubble time is related 
to the age by (cf. equation 5.63 in Peebles (1993). 

If = 0.3 then Hq^ = 0.965 [age]. For — 0.262 the Hubble time would equal the age; in 
other words, the current expansion rate would equal the historic mean. 
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Fig. 1.— 



The appearance of the PixeLens user interface, with all the areas marked. 
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Fig. 2. — Image positions and pixelation for the three lens systems: 1115+080 (upper), 1608+656 
(middle), and 1933+503 (lower). 
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Fig. 3. — Illustration of the Monte-Carlo procedure. The thick lines represent inequality constraints, 
and any point in the pentagonal region enclosed by these lines is an allowed solution. Suppose we 
are at solution point A. We now do the following. (1) Choose a random vertex P of the allowed 
region. (2) Draw the line PQ by drawing PA and extending it to end of the allowed region. 
(3) Choose a random point B on PQ, and move to it. The roles of A and B in this sequence of 
operations is symmetric. Hence the probability of moving move A B equals the probability of 
move B ^ A. 




Fig. 4. — Ensemble-average mass maps for the three lenses. The contours are in steps of k = |, |, . . . 
The large filled circles are the image positions. The small filled circles near the centres are the 
inferred source positions. 
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Fig. 5. — Ensemble-average lens potential for 1115+080. Upper panel: Contours of constant tp, the 
contours shown being the ones passing through the images. Lower panel: Like the upper panel, 
but with the external shear exaggerated by a factor of 3. 
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Fig. 6. — Ensemble-average arrival times for the three image systems in 1933+503: core quad 
(upper), lobe quad (middle), and lobe double (lower). Only the saddle-point contours are shown. 
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Fig. 7. — Arrival-time surface for 1115+080, with contours spaced by 0.005 arcsec^. This represents 
a simple model of the Einstein ring formed by the host galaxy. 
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Fig. 8. — Histogram of the Hubble time from 100 three-lens models. 
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Fig. 9. — Timc-dclay predictions for the core quad in 1933+503: between images 1 and 2 (upper 
panel), between 2 and 3 (middle), and between 3 and 4 (lower). See Saha & Williams (2003) for 
an explanation of how to number images based on arrival times. 
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Fig. 10. — Steepness index a against g, for 1115+080 (upper panel), 1608+656 (middle), and 
1933+503 (lower). The upper and middle panels are analogous to Figures 13a and 16a respectively 
in Williams &; Saha (2000), except for an obvious change of sign. 
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Fig. 11. — Annular density (K)ann against g~^, again for 1115+080 (upper panel), 1608+656 (mid- 
dle), and 1933+503 (lower). 
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Fig. 12. — Mass maps as in Figure 4, but for 1520+530 (upper left), 1600+434 (upper right), 
1830-211 (lower left), and 2149-275 (lower right). 
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Fig. 13.— Histogram of the Hubble time from 100 four-lens models of 1520+530, 1600+434, 1830- 
211, and 2149-275. 
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Fig. 14. — Annular density (K)ann as in Figure 11, but for 1520+530 (upper left), 1600+434 (upper 
right), 1830-211 (lower left), and 2149-275 (lower right). 
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Fig. 15.— 



Schematic representation of how the code in PixeLens is organized. 



